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In order to illuminate the properties of current fluctuations in more than one dimension, we use a 
lattice-based Markov process driven into a non-equilibrium steady state. Specifically, we perform a 
detailed study of the particle current fluctuations in a two-dimensional zero-range process with open 
boundary conditions and probe the influence of the underlying geometry by comparing results from a 
square and a triangular lattice. Moreover, we examine the structure of local currents corresponding 
to a given global current fluctuation and comment on the role of spatial inhomogeneities for the 
discrepancies observed in testing some recent fluctuation symmetries. 
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I. INTRODUCTION 

The understanding of non-equilibrium physics is of 
great relevance for many scenarios ranging from granular 
materials, chemical reactions and molecular motors to 
traffic jams [1-4]. In particular, much interest is directed 
towards the study of the probability of rare events or 
trajectories in stochastic models. This has led to the 
establishment of fluctuation theorems, which are some of 
the most general results for systems out of equilibrium 
(for reviews see [5-7]). Within the stochastic framework, 
interacting particle systems have enjoyed widespread use 
to model non-equilibrium steady states (NESSs). Most 
such models are one-dimensional; we expect a richer phe¬ 
nomenology in more than one dimension just as higher 
dimensional equilibrium systems are qualitatively different 
from their one-dimensional counterparts. 

The study of NESSs in more than one dimension has led 
to the recent discovery of symmetries for global current 
fluctuations in macroscopic systems [8, 9]. In particular, 
by considering a diffusive lattice gas on a d-dimensional 
(hyper-cubic) lattice of side length L, fluctuation relations 
were obtained for a time-averaged global current defined 
as 


J =l [ dr f drj{r,T). (I) 

* Jo Jn 

Here, the local current, is assumed to obey the 

continuity equation, and a diffusive scaling is applied. 
Namely, space is scaled to H = [0,1]^^, and time is scaled 
by 1/L^. Then, according to the macroscopic fluctuation 
theory (MET) the probability of observing a rare global 
current can be calculated from the minimization of an 
action functional that depends on the local value of the 
current and density [10, 11]. Physically, this means that 
out of all the possible ways to generate a fluctuation, the 
overwhelmingly most likely to be realized corresponds 
to a specific optimal density profile (ODP) and optimal 
current profile (OCP). Under some hypotheses, notably a 
spatially homogeneous OCP, it is possible to manipulate 


the MET action functional to obtain that global current 
fluctuations satisfy the relation 


1 , P(J',t) ^ 

hm-j In ———- = E ■ 
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(J-J') 


( 2 ) 


for isometric fluctuations such that 


I J| = \J'\. (3) 

Here, P{J,t) is the probability of observing the global 
current fluctuation J during a time interval t, E is a, 
constant that depends on the boundary and the bulk 
driving of the system, and j • j denotes the modulus. Note 
that this isometric fluctuation relation (IFR) reduces to 
the renowned Gallavotti-Cohen-type fluctuation symme¬ 
try [12-15] for anti-parallel currents but it also relates in 
a surprisingly simple manner currents in different direc¬ 
tions [8]. Furthermore, in [9] the IFR was generalized to 
anisotropic systems, where some discrepancies were also 
noted between the global current fluctuations predicted 
to satisfy the symmetry at a macroscopic level and those 
in models on (large) finite size lattices. Remarkably, these 
fluctuation relations have recently been tested experimen¬ 
tally as reported in [16], where fluctuations of the velocity 
of a tapered rod are shown to be well approximated by 
the anisotropic generalization of the IFR. 

In this paper, we use a continuous-time Markov process 
driven by the boundaries into an NESS in order to study 
the detailed properties of current fluctuations in more 
than one dimension. In particular, a two-dimensional zero- 
range process (ZRP) is solved to study: a) the influence 
of the underlying lattice geometry on the probability of a 
global current (and density) fluctuation, and b) the most 
likely local current structure of the OCP associated to 
a given global current fluctuation. Specifically, we test 
whether the hypothesis of homogeneous OCP may have 
to be adjusted for some systems with open boundary 
conditions, explaining the above-mentioned discrepancies 
observed for rare global currents. 

The paper is structured as follows. In Section H we 
introduce some definitions from large deviation theory 
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commonly used in the study of NESSs. In Section III 
we explain the so-called quantum Hamiltonian formal¬ 
ism which we employ to study the stationary state and 
the probability of measuring rare particle current fluc¬ 
tuations. In Section IV we solve exactly an anisotropic 
two-dimensional (2-(i) ZRP on square and triangular lat¬ 
tices, allowing us to analyse fluctuations and compare 
the effect of the underlying geometry as the system size 
increases. In Section V we refine our calculations to show 
that local current fluctuations have a more complex struc¬ 
ture than implicitly assumed in some other works, and 
highlight the relevance of our results for the anisotropic 
fluctuation relation (AFR) and the original IFR. In Sec¬ 
tion VI we summarize our findings and comment on some 
remaining open questions. In addition, we include various 
technical details in a series of appendices. 

Note that some of these results were already presented 
in a briefer work [9] ; significantly, we here offer convincing 
evidence (Section V) to support the conjecture made there 
on the local structure of current fluctuations. In addi¬ 
tion, we present a new extension to the triangular lattice 
(Section IV) and provide many previously unpublished 
calculational details including a more explicit derivation 
of the AFR than the one shown in [9] (Appendix A) as 
well as the corresponding macroscopic optimization argu¬ 
ment for the 2-d ZRP with non-decreasing interactions 
(Appendix C). 


II. LARGE DEVIATION THEORY: CURRENT 
FLUCTUATIONS 

One of the main goals in this paper is to study the 
structure of the current profiles that yield a particular 
global current fluctuation. However, our results are di¬ 
rectly related to the accuracy of the IFR and the AFR 
for open systems. In this section we introduce these two 
fluctuation relations beginning from the framework of 
large deviation theory. 

The study of NESSs involves understanding the proba¬ 
bility of measuring rare currents. In a lattice-gas model 
for example, a current is understood as the net number 
of particles that jump between two adjacent sites in a 
positive direction (arbitrarily chosen) during a given time 
window. When a system is in an NESS the mean flux 
of particles is, generally, a constant different from zero. 
Moreover, it is known that in many cases such currents 
obey a large deviation principle (LDP), for instance the 
global current J in Eq. (1) satisfies 

e(J) = lim -^lnP(J,t), (4) 

where e( J) is a rate function (RF) encoding the probabil¬ 
ity, P{J,t), of observing a given current in the long-time 
limit [17, 18]. 

In order to calculate the RF, it is useful to compute 


first the scaled cumulant generating function (SCGF) 

e(A) = lim-^ ln(exp (—• J)) (5) 

where A is a vector conjugate to J, and (•) denotes an 
expectation value. It is well known that when the SCGF 
is differentiable we can compute the RF using the Gartner- 
Ellis Theorem, which relates these functions via the Leg¬ 
endre transform [17] 

e(J) = max{e(A) — A • J} . (6) 

As we will see below, much can be learned from the SCGF 
about the probability of the currents, but in the rest of 
this section we remind the reader about some fluctuation 
relations which will be discussed later in the paper. 

At this point it is possible to use the LDP (4) to rewrite 
the fluctuation relation (2) in terms of the RF as 

e(J) - e( J') = E-{J'- J), (7) 

for global currents satisfying |J| = |J'|. Here, the con¬ 
stant E can be seen as an external field driving the system 
out of equilibrium. Moreover, this also implies a symme¬ 
try at the level of the SCGF which is expressed simply 
as 

e(A) = e(A'), (8) 

for values of A such that 

|A-£;| = |a'-f;|. (9) 

Here we note that Eq. (9) also corresponds to the equation 
of a circle for the conjugate parameter of the current, but 
with centre at the constant field E. As mentioned above, 
such a symmetry yields as a special case the Gallavotti- 
Cohen-type relation for forward and backward currents. 
The AFR (derived in Appendix A) is a generalization 
of Eqs. (7)-(9) where ellipses, instead of circles, relate 
current fluctuations in different directions. 

In the following section, we will explain how to study 
fluctuations of a similar space-integrated global current, 
in finite (microscopic) ZRPs. Later, we will explain how 
to scale such a current to compare the results with the 
ones obtained from a macroscopic point of view. 

III. ZERO-RANGE PROCESS: MICROSCOPIC 
APPROACH 

A. Definition of the model 

The ZRP is a model of interacting particles introduced 
in 1970 [19] and since studied on general lattices [20, 21]. 
Particles are allowed to accumulate up to any non-negative 
number on each site of the lattice (e.g.. Fig. 1). The top¬ 
most particle of each site jumps to a neighbouring site 
after an exponentially distributed waiting time, where 
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the hopping rate is proportional to an on-site particle 
interaction factor, Wn- As the name of the model suggests, 
Wn depends exclusively on the total occupation of the 
departure site. Indeed, such an interaction can cause 
a phase transition where a macroscopic proportion of 
particles in the system accumulates on a single site of 
the lattice [22, 23]. Similar condensation phenomena are 
of wide interest in connection with granular systems [24] 
and wealth models [25] among other topics. 

In order to study this zero-range model we employ 
a general framework [26], referred to as the quantum 
Hamiltonian formalism, in which the master equation of 
the system is written in a form resembling a Schrodinger 
equation. Within this approach one can compute the 
probability of particle configurations in the system, as well 
as other important quantities, such as the time-averaged 
particle current. 

We begin by defining the conhguration, n = 
(ni, n 2 ,..., njv), containing the number of particles on 
each of the N sites of the lattice at a given time. Then, 
each configuration n is associated with an element of a 
basis, jn), to construct the probability vector 

lP)=^P(n)ln), (10) 


where we explicitly show the Hamiltonians corresponding 
to the square and triangular lattices (Eqs. (Bl) and (B3) 
respectively). We now turn to study the time-independent 
solution of (11), i.e., the steady state. 

B. Steady state solution 

Typically, to drive a system out of equilibrium, we let 
it interact with more than one reservoir and expect it 
to reach an NESS in the long-time limit. In the present 
context, by considering a ZRP with open boundary con¬ 
ditions where the input and output rates are different 
at the two borders we expect the system to approach a 
time-independent solution with constant mean current 
through the system. This means that the left-hand side 
(LHS) of (11) vanishes, leaving us with the eigenvalue 
equation 

H\P*)=0 (13) 

implying the stationary state, \P*), is the eigenvector 
with eigenvalue zero. Similarly to certain other interacting 
particle models, for the ZRP with open boundaries the 
vector jP*) is given by the product measure [23, 27] 


where P{n) is the probability of finding the system in 
configuration n. 

The time evolution of the probability vector is described 
by the master equation 


d\P} 

dt 


-H\P). 


( 11 ) 


Here the stochastic generator, P, also called the Hamil¬ 
tonian, contains the hopping rates between all states of 
the system, and can be written in terms of the ladder 
operators 


/ 0 0 0 \ 


O 

o 

1 0 0 ... 


0 0 1^2 ■ • ■ 

0 1 0 

and a~ = 

0 0 0 

V ^ 


1 ^ 


( 12 ) 

which act exclusively on the component of the config¬ 
uration vector. 

On each lattice shown in Pig. 1, particles from the bulk 
jump with rates pk and qk in the positive or negative k- 
direction (respectively), as indicated by the insets. Thus, 
a particle jump is represented in the quantum Hamil¬ 
tonian formalism by the simultaneous annihilation and 
creation of one particle, at the departure and target sites 
respectively, with the operators (12) times the correspond¬ 
ing hopping rate. Eurthermore, particles are injected at 
constant rates at and 5k, or extracted with rates jk and 
Pk, both at the left and right boundaries respectively. 

Note that for the two dimensional systems we study, it 
is convenient to identify sites, and corresponding ladder 
operators, with two subindices as done in Appendix B 


|P*) = |P*)®|P*)®...®|P^), (14) 

i.e., the probability distribution factorizes over the sites. 
For the ZRP, it can be shown that the marginal for the 
site is 


(15) 

Tli 

where the probability of finding particles on site i is 
given by [23] 


fli 

P:{n{) = Z-^zT\[w-\ (16) 

Here is the fugacity of site i and Zi is the grand 
canonical partition function 


oo 3 

(ii') 

0 n—1 

Note that for some choices of the interaction term Wn the 
radius of convergence, Zmax, of the sum (17) may be finite. 
Within the range of values where the partition function is 
well defined, the site densities are related to the fugacities 
via the equation 


Outside this range, the diverging Zi corresponds to the 
accumulation (condensation) of infinitely many particles 
on site i. Here we aim to study current fluctuations 
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FIG. 1. (Color online) Different underlying geometries, (a) square lattice and (b) triangular lattice, with periodic boundary 
conditions in the vertical direction (j/-direction) and open boundary conditions in the horizontal direction (^-direction). Sites 
are identified by their row (^-coordinate, first subindex) and column (x-coordinate, second subindex) as shown in the insets. 
In the triangular geometry, sites on the same row are joined by a zigzag line. In the microscopic model open boundaries are 
specified by the input rates (q and S) as well as the output rates (7 and /3) for left and right boundaries respectively (see also 
Fig. 7); the input rates can be related to the left and right reservoir densities pi and pr- 


within the fluid regime of the ZRP (i.e. no condensation), 
for this purpose, it is sufficient to consider Wn as an 
increasing function of the number of particles. As in the 
one-dimensional ZRP with open boundaries [23], it turns 
out that the fugacities are independent of Wn- However, 
Eq. (18) will become relevant when we study the density 
profile related to a given current fluctuation [28]. 

In practice, to compute the fugacities we note that the 
creation and annihilation operators affect only the corre¬ 
sponding site component of the stationary state eigenvec¬ 
tor according to 

at\P:) = z-H,\P:) (19) 

c^\p:) = z^p:). ( 20 ) 

Here, we have defined the diagonal matrices di with ele¬ 
ments djk = WjSj^k where icq = 0 by definition and Sj^k 
is the Kronecker delta. Then for a lattice of N sites, 
Eq. (13) can be reduced using (19) and (20) to a system 
of N coupled equations for the fugacities of the system. 
Note that in this framework, conservation of probability 
implies that the corresponding left eigenvector has every 
component equal to unity; we denote such an eigenvector 

by (i|. 

C. Current fluctuations 

In addition to the probabilities of configurations, it is 
also possible to study the probabilities of fluctuations 


of particle currents within the system. Specifically, our 
goal is to quantify the probability that the time-averaged 
number of particle jumps between nearest neighbouring 
sites, in the whole lattice or a subset of it, attains a 
given rare value. This means that we have to modify the 
quantum Hamiltonian to count the number of particles 
that jump within the lattice during the observation time 
interval. 

To avoid confusion with the macroscopic current, we 
will use the variable I to denote the space-integrated 
microscopic current measured across the lattice. Later 
we will clarify how to rescale this current to compare 
it with the macroscopic approach, but first we explain 
how the current is constructed. We define a certain 
time evolution of the system as the set of configurations, 
{cr} = {cti, ( 72 , ■■■, O',-}, visited by the system during the 
time interval [0,t]. In one dimension it is clear that the 
net number of particle jumps is counted with an antisym¬ 
metric function; we let Oai+i,ai take the value - 1-1 when 
particles jump forwards and —1 when particles jump back¬ 
wards anywhere in the lattice (for a more general case see 
e.g. [29]). This way, the space-integrated current in one 
dimension is defined as the sum 

1 

1'^}) ~ Qg-i+i ,o-i ■ (21) 

i=0 

In higher dimensional lattices, we will be interested in 
a similar vectorial variable with component Ik(t, {cr}) to 
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count the number of jumps in the /c-direction. For now we 
keep the one dimensional notation in order to demonstrate 
the framework. 

In analogy to the macroscopic case, to compute the 
microscopic RF 

^l{I) = lim -ilnP(/,t), (22) 

>-oo t 

it is useful to define first the SCGF 

s-lM = lim — - ln(exp (—<A/)). (23) 

t—foo t 

It can be shown that the average on the right-hand side 
(RHS) of this relation can be written as {exp{—tH)), 
where iJ is a modified Hamiltonian. In order to obtain 
H, we have to multiply the terms of H corresponding 
to particle transitions by exp(—A) for jumps in the pos¬ 
itive direction and by exp(A) for jumps in the negative 
direction [30]. 

In cases when the eigenvalue spectrum of H is gapped 
the calculation of the SCGF can be done by noticing 
that in the long-time limit the exponential of the lowest 
eigenvalue, dominates the average (23). This leads 
to the result 


eL(A) = C(A), (24) 


where we have assumed that the prefactors arising from 
the eigenstate decomposition are finite. Breakdown of 
this condition signifies condensation. 

Furthermore, the right ground-state eigenvector, |'0), 
turns out to have the same form as the stationary 
state (I4)-(17) but with some modified fugacities, Zi(A). 
In principle it is possible to calculate exactly the modified 
fugacities by using relations analogous to Eq. (19) 
and Eq. (20), allowing us to determine also the SCGF. No¬ 
tice that the lowest eigenvalue does not vanish in general. 
Indeed, one can verify consistency with the stationary 
state by substituting A = 0, for which the eigenvalue does 
become zero. 

Finally, the RF is calculated via a Legendre transform 
similar to Eq. (6) for microscopic currents. We remark 
that when the transform cannot be computed analytically, 
we can use the implicit relations 


deL{\)) 

dX 

deL{I) 


(25) 


to calculate it numerically. 

To obtain the density profile which gives rise to these 
currents, we have to compute the mean local occupation 
{rii) taking into account the dynamics of the modified 
Hamiltonian. To do this, we need both left and right 
modified eigenvectors corresponding to the ground state 
eigenvalue. The left (row) eigenvector, {ipl, is again a 
product with terms denoted by (V^ij = (I, zf,...). To 


calculate the left fugacities, Zi(A), we use the modified 
Hamiltonian on ('0|, where the ladder operators act ac¬ 
cording to the relations 

{ipi\af = {i^i\zi (26) 

{A\a~ = {'ipilz-^di. (27) 

Here the dependence on A is left implicit; this is done 
from now on for both Zi and z^. 

Using Eqs. (26) and (27), the left fugacities are ob¬ 
tained in terms of the lattice parameters as outlined in 
the framework above for the right eigenvector. Note that 
the components of the left and right eigenvectors are not 
the same in general. Moreover, the typical density at site 
i associated to a given current fluctuation can now be ex¬ 
plicitly calculated using the definition Pi{X) = {'ipi\ni\'ipi), 
where rii is the diagonal operator for the number of parti¬ 
cles on site i. This leads to a relation between densities 
and fugacities similar to the grand canonical identity (18) 
for the steady state, but with the replacement of Zi by 
ZiZi. In particular, for the interaction Wn = n the site 
density is determined by Pi(A) = ZiZi which reduces in 
the stationary state to / 9 i( 0 ) = Zi since Zi equals unity for 
A = 0. For bounded the product (V’1'0) might diverge 
which again generically indicates condensation [30]. 

In the following we show how to extend the formalism 
presented here to study the ZRP on different 2-d lattice 
geometries. In particular, we need to study the influence 
of the underlying lattice structure on current fluctuations 
in large finite lattices. To do this, we will modify the 
Hamiltonians with factors exp(=pAfe) to count simulta¬ 
neously the number of jumps along the corresponding 
positive or negative /c-directions. 

Specifically, the modified Hamiltonians for the systems 
shown in Fig. 1 are given explicitly by Eq. (B2) (with 
Xy — 0) for the square lattice, and by Eq. (B4) for the 
triangular lattice. In this manner, for fixed lattice sizes, 
we will compute the SCGF 

es'(A) = lim — - ln(exp (—tA •/)), (28) 

t —>00 t 

for the square geometry and 

eT(Ai, A 2 , A 3 ) = lim ln(exp (-t(Ai, A 2 , A 3 ) 

t (29) 

• {Ii, 12,13))), 

for the triangular geometry. To avoid confusion, from 
now on we will use bold characters to denote vectors in 
Cartesian coordinates, and we will write explicitly the 
components of the variables in the triangular lattice. 


IV. EFFECT OF LATTICE GEOMETRY ON 
CURRENT FLUCTUATIONS 

It is known that in some cases the underlying geometry 
in lattice-based models can have a significant effect on the 
results obtained from a microscopic point of view. For 
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example in a different context, at the level of the univer¬ 
sality of phase transitions in equilibrium systems, some 
differences were found numerically between square and tri¬ 
angular lattices [31] and later predicted theoretically [32]. 
Here, we use the quantum Hamiltonian formalism to cal¬ 
culate the RF of the global current in the two geometries 
shown in Fig. 1 and thus investigate the influence of the 
lattice on the dynamical properties of the ZRP. 

One reason we are interested in the RF of particle cur¬ 
rent fluctuations and their associated ODP is to confirm 
that the expressions from both lattices converge to the 
same function and recover the hydrodynamic result under 
the appropriate scaling. We shall give more details of 
the scaling in subsection IV C, but first we calculate the 
SCGF for finite lattices using the microscopic approach 
introduced above. 


and site exit rates 


3 


Ri = ^^Vk + + 72 + 73 

k^l 

3 

(37) 

Rr — qk Y P\ P 2 

k^l 

3 

(38) 

R = ^{PkY qk) ■ 

(39) 


It can be checked that the same difference equation and 
boundary conditions are obtained for the square lattice 
with coefficients: 

P = Rr = P + Qx + Py + Qy 

Q = q.e^- Y = pye-^y+qye>'y 

Tt — Px “t“ qx “t“ Py qy 


A. General solution on square and triangular 
lattices 

As mentioned above, to calculate the probability of 
fluctuations of the global current, we have to measure 
the number of jumps throughout the lattice in the time 
interval [0,t]. The modified eigenvector, \'ip), associated 
to the lowest eigenvalue of the modified Hamiltonian, 
H, obeys relations analogous to (19) and (20) with the 
modified fugacities z^i. 

For each lattice we apply the corresponding H, (B2) 
or (B4), to and use the eigenvector condition to obtain 
that the coefficients of the matrices dj^i in the resulting 
expression have to vanish. This leads to the recursion 
relation for the modified fugacities of the right eigenvector 

Q^j,i+i d- {Y — R)zj^i + Pzj^i-i = 0 (30) 

with boundary conditions 

Qzj,2 + {Y - Ri) Zjp + Ai={) 

Ar + (Y — Rr) ZjX + PZjX-1 = 0, 

for the left- and right-hand side respectively. Here, the 
uppercase parameters for the triangular lattice correspond 
to the effective bulk rates 


P = P2e ^^+P3e (32) 

Q = 926 ^" + 936 ^" (33) 

r =pie-^i+gie^b (34) 

boundary injection rates 

2 

= (35) 

k^l 

2 

A. = (36) 

k^l 


Rl =Px+ 7 +Py+qy Ar = . 

Here we have omitted the subindices of the boundary 
parameters as particle jumps in and out of the system 
contribute only to the current in the a;-direction. 

We make use of the periodic boundary conditions to 
argue that we only have to solve Eq. (30) and (31) for 
a single row. In other words, since the fugacities are 
invariant in the y-direction the system can be treated as a 
quasi-one-dimensional chain. The difference equation (30) 
can be solved exactly, but the expressions are too cum¬ 
bersome to handle; in practice we use a computer algebra 
package to calculate exact numerical values for the mod¬ 
ified fugacities. An analogous calculation is required to 
find the components of the left eigenvector, (^], using 
relations (26) and (27). As we shall demonstrate, the 
seemingly technical analysis of this eigenproblem allows 
us to investigate both the probability of given current 
fluctuations (since the eigenvalue generically gives the 
SCGF) and the mechanisms leading to them (since the 
modified fugacities can be related to densities). 

The results of the following subsections are based on 
the fact that the ground state of the modified Hamiltoni¬ 
ans (B2) and (B4) for these lattices can straightforwardly 
be written in terms of the modified fugacities such that 
the SCGF for the triangular lattice is given by 

3 

er(Ai,A2,A3) = l'^ {au +4 - Pke~^’‘ZL) , 

fe=2 

(41) 

whereas for the SCGF of the square lattice we have 

es{X) = L [a + 5 - - pe~^^ZL) ■ (42) 

Note that due to the symmetry imposed by the periodic 
boundary conditions in the y-direction only the second 
subindex, related to the ^-direction, is needed to identify 
the fugacities. Finally, as mentioned above, we can calcu¬ 
late the average (V'lnij^) with the modified eigenvectors, 
which corresponds to the ODP from a microscopic point 
of view. 
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B. Matching diffusive processes in square and 
triangular lattices 


geometries can be obtained by noticing that currents in 
the triangular lattice have components 


Before we can compare the solutions obtained from the 
square and triangular lattices, it is necessary to choose 
carefully the bulk and boundary hopping rates in order 
to achieve an equivalent behaviour in both systems. 

Firstly, since we are interested in modelling diffusive 
dynamics in the hydrodynamic limit, we consider sym¬ 
metric hopping rates qk = Pk- Additionally, we match 
the extraction boundary rates to the bulk hopping rates 
Ik = Pk = Pk SO that the boundaries act simply as reser¬ 
voirs. 

Now to obtain a mapping for the bulk hopping rates 
between the triangular and square lattices, we need to 
equate the particle transport bearing in mind the lattice 
spacing. In other words for our choice of diffusive dy¬ 
namics, we have to equate the mean square displacement 
in both lattices. Mathematically, this implies that the 
hopping rates of the triangular lattice are mapped to the 
square lattice via the relations 


jx = j 2 cos p + ja cos (j) 

3y = ji + h sin p - js sin p. 

Then, we can use the chain rule on Eq. (25) together with 
the relations (46) to obtain 

Al = Xy 

X 2 = Xx cos <j) + Xy sin (j) (47) 

A 3 = Xx cos 4> — Xy sin 0 , 

which are the appropriate conjugate variables to compare 
the number of particle jumps on a triangular lattice with 
those in a square lattice. 


C. Hydrodynamic limit and optimal density 
profiles 


Px = P2 COS^ p + Ps COS^ (j) 

Py = Pi + P2 sin^ 4> + P3 sin^ 


(43) 


where p = tt/G. 

To obtain diagonal matrices for the diffusion and mo¬ 
bility coefficients as required for the process on the square 
lattice, we need to identify ps = P 2 - Such a choice of 
hopping rates leads to the simplified mapping 


Pi =Py- 
‘3‘Px 

Y 


Px 

3 


or 


P2 = 


3p2 

Pv=Pi + -J- 


(44) 


Now we can check that for an isotropic choice of hopping 
rates in the square lattice (i.e., Py = Px), our transforma¬ 
tions yield isotropic rates in the triangular geometry. Simi¬ 
larly, we can use (44) to confirm 2{px+Py) = 2{pi+p2+P3), 
so the exit rate from a bulk site is the same in both lat¬ 
tices. 

The same reasoning can be used to determine the map¬ 
ping of the boundary rates, which yields the analogous 
expressions 


2ci 25 i AK\ 

^k — 5k — (^^) 

These relations conserve the injection-extraction ratios, 
Ci/Px = OL 2 IP 2 = OLzjps (and analogously for the RHS 
boundary), which is equivalent to reservoirs with the same 
fugacity, zi = ak/pk and Zr = 5k/pk, for both lattices. 

Using the transformation relations (44) and (45), parti¬ 
cle diffusion on the two lattices can be related. To convert 
current fluctuations in the triangular lattice to Cartesian 
coordinates, we have to specify how to count particle 
jumps with the quantum Hamiltonian formalism. An 
appropriate relation between the triangular and square 


In this subsection, we focus on explaining the scaling of 
our results obtained from the microscopic approach with 
the goal of determining the influence of the underlying 
lattice geometry for large systems. We will also compare 
the microscopic results with those obtained in Appendix C 
using the MET. We begin by discussing the scaling for 
the SCGF (28) of the square lattice, as it is more intuitive 
than the triangular geometry which will be explained 
immediately after. Later, we will explain how to scale the 
density profiles, which is done in a similar way. 

As hinted in the introduction, to obtain the hydro- 
dynamic limit of diffusive systems starting from a mi¬ 
croscopic approach, a spatial and temporal rescaling is 
needed. Specifically, space is scaled as I/L and time as 
1/L^, where L is the linear size of the microscopic system 
(see e.g. [33]). In the quantum Hamiltonian formalism, 
this leads to dividing the conjugate parameters by the 
number of bonds in the corresponding direction, as we 
measure particle jumps throughout the lattice. In d di¬ 
mensions, the temporal and spatial rescaling also combine 
to give a factor of L^~‘^ multiplying the SCGF [34] but 
this reduces to unity for our square 2-d system. Hence 
we expect the macroscopic SCGF, e(A), to be given by 
the limit 


e(A) 


lim 

L—>oo 


L + l 


L 


es 


Xx Xy \ 

l + vy) ■ 


(48) 


Here we have included an additional factor of (L -I- 1)/L 
to remove finite size effects in the cc-direction of small 
lattices; the large-L limit is clearly unaffected by this. 

To obtain the correct scaling for the triangular lattice, 
we have to remember that the length of the lattice in the 
cc-direction is smaller than in the square lattice, being 
multiplied by a factor of cos (f). This can be compensated 
by a modified-length spatial and temporal scaling, leading 




US to the limit 


0.15 


e(A) 


L 1 / Ai 

lim 7-7- er j - 

L-s-oo L COS^ (p \L cos (p 

A2 A3 

[L + 1) cos 0 ’ (L + 1) cos (j) 


(49) 


where Ai, A 2 and A 3 are given by (47). We have checked 
that the additional cos (j) factors in the argument of ex can 
be removed by considering a lattice with Lj cos(f> sites in 
the a:-direction (i.e., spatial length L) and L sites in the y- 
direction. However, this produces more complicated finite 
size effects since the number of sites has to be rounded 
to an integer. 

In the present work, we show the results for the inter¬ 
action Wn = n, but we have checked also the case with 
Wn = w (w constant) within the fluid regime which leads 
to similar findings. Indeed, the SCGF is invariant with 
respect to the interaction as long as there is no conden¬ 
sation [35], but the relation between the density and the 
fugacity (and hence between boundary rates and reser¬ 
voir densities) does change. The special case of = n 
is particularly illuminating because densities then turn 
out to be proportional to fugacities, so calculations of 
the latter offer direct physical insight into the optimal 
profiles. In addition to verifying the hydrodynamic limit 
with the above scaling, studying the SCGF also provides 
a convenient way to test the AFR from a microscopic 
point of view. 

In Fig. 2 we plot the RHS of (48) and (49) for increasing 
lattice sizes and values of A for which Eq. (A16) is satisfied. 
We assume bulk and boundary hopping rates. 


a = 1/2 
7 = /3 = 1 
5 = 1/10 


Vx — Qx — 1 

Py ~ ~ 


(50) 


which make particle diffusion anisotropic; of course, it 
is also possible to test the IFR for isotropic rates. The 
rescaled microscopic results are compared with the numer¬ 
ical Legendre transform of the macroscopic RF obtained 
in Appendix C with reservoir densities pi = a and Pr = S. 
We can see that both microscopic SCGFs converge to the 
same function when L —>■ 00 . As might be expected (due 
to a larger number of bonds), with the triangular lattice 
the SCGF has a quicker convergence towards the hydro- 
dynamic limit than with the square lattice. However, it 
can also be observed that this limit does not agree with 
the result obtained using the MFT under the assumption 
of homogeneous OCPs. 

Turning our attention to the AFR, note that in Fig. 2 
we parametrize in polar coordinates (with angle 9) the 
values A for ellipses centred at 


E = 


1 

2 




(51) 


on which the SCGF is predicted to be constant. As shown 
in [9] , here the AFR is satisfied by the macroscopic results 
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FIG. 2. (Color online) SCGF for \{0) on concentric ellipses 
around E = 1/2 (ln(a/3/(75)), 0) with principal axes in the 
a;-direction of length = {0,0.266,0.533,0.8} from top to bot¬ 
tom. Square lattice with L = {6 (a), 10 (•), 10® (•)}, triangu¬ 
lar lattice with L = {6 (A), 10 (O), 10® (□)}, and macroscopic 
approach (solid line). Hopping rates given by Eq. (50). 


for all points on the ellipse but only for certain angles from 
the microscopic point of view. Importantly, agreement 
in the hydrodynamic limit of the microscopic SCGFs 
of the two lattices, indicates that the discrepancies are 
not related to the underlying structure. We continue to 
investigate further the ODP from both approaches. 

The microscopic ODP is obtained as explained 
in Section HI, but also have to be rescaled before 
we can compare them with the macroscopic ODP 
computed in Appendix C. Firstly, notice that within 
the modified Hamiltonian dynamics the ODP at site 
i is given by the mean occupation, (n^), and has to 
be compared with the macroscopic ODP at x = ijL. 
Secondly, scaling of the conjugate variables is done 
the same as in the arguments of (48) and (49). This 
way, we are able to compare ps {Xx/{L + 1), Xy/L) and 
Pt (Ai/(Lcos (/), A 2 /((L -I- l)cos(/), A 3 /((L -I-1) cos </>)) 
for the square and triangular lattices, as well as p ( J(A)) 
for the macroscopic approach. 

Since the ODP has no dependence in the y-direction 
(due to periodic boundary conditions), in Fig. 3 we plot 
the x-projection of the profiles for two values of A which 
are predicted to satisfy the AFR. We have used the same 
bulk and boundary parameters as for the SCGFs plotted 
in Fig. 2. It can be seen from the upper panel of Fig. 3 
that our calculations for both lattices match closely the 
macroscopic solution (solid line) for A in the x-direction 
even for lattices with L = 10. In the lower panel we 
choose a current in the y-direction only, specifically a 
current of appropriate magnitude such that the macro¬ 
scopic ODP remains invariant with respect to the upper 
panel. The microscopic-approach ODPs from the two 
lattices converge towards the same function in the hy¬ 
drodynamic limit, but not to the MFT prediction. This 
again indicates that the AFR is not exact between these 
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X 

(a) 



X 

(b) 


FIG. 3. (Color online) Optimal density profiles on L x L-site 
lattices. Hopping rates: a = 1/2, <5 = 1/10, Px = 1, Py = 1/2, 
and Wn = n. ODP from microscopic approach on square lattice 
(□), and triangular lattice (A) with symbols of decreasing 
size for L = {10,20,10®}; ODP from macroscopic theory 
with blue solid line, (a) Current fluctuation in i-direction: 
A ~ (—0.6953,0), i.e., J ~ (0.9523,0). (b) Current fluctuation 
in i/-direction: A ~ (0.8047,-2.1213), i.e., J ~ (0,0.8928). 


current fluctuations. We suggest that the assumption 
of a space-homogeneous OCP in the MFT lies behind 
this discrepancy, and we will investigate it further in the 
following section by looking for a more detailed structure 
of the current fluctuations in L x L square lattices. 


V. STRUCTURE OF OPTIMAL CURRENT 
PROFILES 

In this section, we extend our study of global current 
fluctuations, to gain a fine-grained understanding of the 
underlying local structure. Specifically, we seek informa¬ 
tion about the OCP giving rise to a particular global 
current fluctuation. In contrast to hypothesis in) of the 
APR (see Appendix A), we anticipate finding some spa- 



FIG. 4. (Color online) We count particle jumps in the y- 
direction within the highlighted region of width e and calculate 
the joint SCGF e{Xy,\). 

tial dependence (with a similar structure expected for 
isotropic systems as relevant for the IFR). This conjec¬ 
ture can be motivated by remembering the definition of 
the global current and the implications of measuring a 
rare realization of it. When we calculate the RF of a cer¬ 
tain fluctuation J, what we are considering is a space- and 
time-average of the number of particle jumps throughout 
the lattice. However, there could be many local current 
profiles, with different spatial dependence, leading to this 
average. From all such profiles we want to find the OCP. 

In order to gain a deeper understanding of the fluctua¬ 
tions in the 2-d ZRP on a square lattice, we consider the 
joint probability distribution function (PDF) of a global 
current and a local current in the y-direction of a vertical 
strip, V, as indicated in Fig. 4. The relative area of V is 
kept constant for all lattice sizes. This implies that the 
width of V is made proportional to the lattice length as 
we increase the number of sites in the system. Thus we 
anticipate that properly rescaled microscopic results will 
approach a consistent hydrodynamic limit for increasing 
L. For the purposes of discussing this limit we use the 
macroscopic notation, but the underlying calculations are 
still done using the microscopic approach. 

To look for the structure of the OCP, we compute the 
joint RF of the local and global currents, 

e{Jy{xo),J) = . (52) 

Here, J is defined as in the MFT according to Eq. (1) 
and the local current corresponds to 

Jy{xo) = \ [ dr [ drjy{r,T) (53) 

c Jo Jv 

where Xq is the left boundary of V. We keep a fixed value 
of J and move the location of V along the x-direction 
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in order to capture the statistical behaviour of the local 
current Jy{xo) in a more detailed way. 

As before we first calculate the microscopic SCGF as 
the lowest eigenvalue of a modified Hamiltonian; here 
counting the number of particle jumps in V along the y- 
direction requires us to modify some terms of the stochas¬ 
tic generator with the additional variable Ay, where the 
dependence of this modification on xq is left implicit 
throughout. Specifically, calculating the modified fugac- 
ities involves solving a recursion relation similar to (30) 
and (31), where introducing Ay imposes the new relation 

Qzj^i+i {Y — R)zj^i -f Pzj^i-i = 0, (54) 


for sites within V. Here Y = . 

Taking into account this modification, we were not able 
to find an exact analytical expression for the SCGF, but 
we still obtain a complete system of linear equations which 
can be solved numerically. Finally, the SCGF is rescaled 
similarly to (48) to obtain the hydrodynamic limit. 

To probe the spatial dependence of the current profile. 


we begin by considering how the RF, ej yJy{xo)j ■= 

e^Jy(xo),J^, with fixed global current J = {Jx,0), 

changes as V sweeps the lattice. This is equivalent to calcu¬ 
lating the RF of the conditional probability of measuring a 
local current, given a fixed global current weighted by the 
probability of that global current, (i.e. P{Jy{xQ)\J)P{J))\ 
the RFs of the conditional and the joint probabilities dif¬ 
fer only by a term independent of Jy{xo). In practice, we 
will focus on studying the corresponding SCGFs. 

In Fig. 5 we plot the joint SCGF of Ay and fixed A, 
ex(Xy). From the top panel (where A is fixed in the x- 
direction) we can see that at the three chosen positions 
on the lattice, we have dex{Xy)/dXy\j^ = 0, implying 
that the local mean current in the y-direction vanishes 
in all cases. Additionally, the SCGF becomes broader as 
V approaches the right boundary (broader SCGF means 
that the absolute value of the second derivative is smaller). 
Taking into account that the variance of the local current 
can be calculated as —d‘^ex{Xy)/dXy\i^^^, implies that 

Jy{xo) is less prone to fluctuations near the right reservoir. 
This can be understood physically as having a higher 
chance to see variations of the current where the sites 
have more particles available, as long as the system is in 
the fluid state. The fact that the variance of Jy{xo) is 
spatially dependent, means for its conditional PDF that 
in general when xq a;g 



(a) 



(b) 


FIG. 5. (Color online) ex(Xy) for A’s predicted to sat¬ 
isfy the AFR, and local conjugate parameter. Ay, for a slit 
of relative width e = 1/20 with its left boundary located 
at xo = {0,2/5,4/5} (parabolas of increasing broadness); 
rescaled microscopic results from a lattice with L = 10®. Same 
boundary rates, bulk hopping rates, and as Fig. 3. (a) Cur¬ 
rent fluctuation in ^-direction: = 0.1 and 6 = n, i.e., 

A ~ (0.7047,0). (b) Current fluctuation in diagonal direction: 

= 0.1 and 6 = 57r/4, i.e., A ~ (0.7340, -0.1). 


maximum of the SGGF displaced. One can easily check 
that the displacement corresponds to 

Ey = —Ay, (56) 

which can be seen as an artificial field in the driven dy¬ 
namics caused by the global conditioning. Significantly, 
taking the derivative of the SCGF we have 


P(Jy(xo)|J)^P(Jy(x'o)|J) (55) 

even if J = {Jx,0). 

In the same manner, we calculate the SCGF assuming 
a fixed global current fluctuation away from the x-axis 
(9 = 57r/4). The result is shown in the bottom panel 
of Fig. 5. In this case, we recover the same behaviour 
as before for the broadness of the SCGF, but with the 


E[Jy{Xo)\J[ 


dexCXy) 

dXy 



(57) 


which is no longer constant at different locations of V, 
meaning that in general when xq ^ Xq 


E[JyiX0)\J]^E[Jyix'0)\J]. 


(58) 
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Here, E['|-] denotes the conditional expectation of a local 
current. The inequality (58) implies physically that for 
a specific global current, the average local current at 
xo and at Xq is not the same, causing the OCP of the 
corresponding J to be inhomogeneous. 

To see in more detail the implication of the inequal¬ 
ity (58), we plot in Fig. 6 the local mean current in the y- 
direction for different global current fluctuations predicted 
to satisfy the AFR (A13). Namely, we fix A on ellipses 
centred at the constant field (51) and obtain E[Jy(xo)| J] 
along the lattice. In particular, we have chosen A at an¬ 
gles 9 = {tt, 77 r/6, 57r/4,47r/3,37r/2} belonging to ellipses 
where the distance to the centre in the x-direction (i.e., 
at 0 = 0) is Tj, = {0.1,1.5}. By taking values in the lower 
left quarter of the ellipse around E, we obtain currents in 
the upper right plane after being mapped by the Legendre 
transform (as reported in Fig. 6). Here, we can see that 
the mean current is homogeneous only when a fluctuation 
of the global current is precisely in the x-direction. 

Indeed, the inhomogeneity of the local current is con¬ 
sistent with the action functional of the macroscopic fluc¬ 
tuation theory (see Eqs. (A2) and (A3)): notice that 
the RF is inversely dependent on the mobility coefficient, 
cr(p). In particular, in the ZRP with interaction = n, 
<j{p) oc p so the mobility is an increasing function of the 
density which implies that it is more cost-effective for 
the system to generate a current fluctuation where it has 
a higher density, typically near the left reservoir. This 
argument can also be made for other ZRPs where the 
mobility coefficient increases with the density, e.g., with 
interaction Wn = w (constant w) the mobility coefficient 
is related to the density according to a{p) oc wp/{p + 1). 

We have also checked that adding the mean value of 
the local currents on disjoint regions V (covering the 
whole lattice), is consistent with the value fixed for the 
y-component of the global current. This indicates that by 
considering smaller widths the local mean current profile 
should converge to the OCP. The analysis of this section 
therefore suggests that the OCP can be space dependent, 
which could be responsible for the discrepancies between 
the microscopic and (current-homogeneous) macroscopic 
approaches, as well as the fact that the AFR (and the 
IFR) is not exactly satisfied for currents in the j/-direction. 


7.x 10"^ 

6.x 10“^ 

5.x 10"^ 

"A 4.x 10"^ 
iK* 3.x 10"^ 

2.x 10"^ 
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0 . 
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FIG. 6. (Color online) Local mean currents in y-direction 
as a function of xo for fixed J predicted to satisfy the 
AFR. Symbols show numerical values for V of relative width 
e = 1/100 at xo G {0,1/5,2/5,3/5,4/5,99/100} rescaled 
from a lattice with L = 10^. Solid lines show interpola¬ 
tion with fourth degree polynomial. Boundary and bulk 
hopping rates given by (50). Global currents at angles 
6^ = {0 (■), 7r/6 (a), 7r/4 (#), tt/S (T), 7r/2 (•)} on ellipses sat¬ 
isfying (A14) passing through (a) J ~ (0.0448,0) and (b) 
J ~ (0.9522,0). 

VI. DISCUSSION AND OUTLOOK 

We have studied the ZRP on square and triangular lat¬ 
tice geometries calculating exactly the fugacities through¬ 
out the lattice, the SCGF for global current fluctuations, 
and the density profiles associated to such fluctuations. 
We also used these results to test a recently predicted 
symmetry for anisotropic systems (the AFR). Since the 
ZRP we studied is solved analytically, our results have an 
advantage compared to other studies of the same class of 
models where numerical simulations are needed to test 
convergence towards macroscopic predictions. For ex¬ 
ample, in [8] the IFR (for isotropic systems) was tested 
using the Kipnis-Marchioro-Presutti (KMP) process and a 
hard-disk fluid. In particular, despite using an efficient al¬ 
gorithm, the KMP process was simulated for a maximum 
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lattice size of L = 32 (i.e. 32^ sites). 

In [9] we raised as an observation that given the large 
lattice sizes considered, up to L = 10^, the SCGF and 
ODPs obtained from the microscopic approach did not 
seem to converge exactly to the macroscopic prediction. A 
similar result was obtained here for the triangular lattice, 
with a quicker convergence of the microscopic results 
to the same hydrodynamic limit as observed with the 
square geometry. We believe the discrepancy with the 
macroscopic prediction is caused by the fact that, in the 
MFT, the OOP was assumed to be spatially homogeneous. 
Our analysis here of the local structure of the current 
fluctuations (within the quantum Hamiltonian formalism) 
indicates that such an invariance does not hold in general. 
In fact, one could also relax the assumption of space 
homogeneity in the MFT and it would be interesting 
to check the resulting OOP, ODP, and SCGF with the 
hydrodynamic limit of our results. (Significant work in 
this direction has already appeared since the submission 
of our original preprint [36].) 

Crucially, the hypothesis of spatial invariance of the 
OCP is also used in the derivation of the AFR (and the 
original IFR). The finding of spatial inhomogeneity thus 
explains the fact, that for fluctuations away from the field 
direction, the AFR does not hold exactly in this model. 
However, we still expect some kind of fluctuation relation 
along the lines of the AFR without assuming homogeneous 
OCPs. Such a generalization would presumably not have 
the same simple structure of (A13) and (A14) but relate 
only local rotations of the current (compare with the 
discussion for the IFR in [8, 34]). We emphasize that the 
usual AFR is still significant for experiments [16], because 
it is a good approximation for fluctuations close to the 
forward direction and therefore enables the testing of 
fluctuation symmetries without the need to measure rare 
backward fluctuations. Furthermore, for systems with 
periodic boundary conditions in every direction, the OCP 
is not expected to have any local structure, so this type of 
spatial fluctuation relation should be exactly satisfied [37] . 

Finally, we point out that knowing the local structure 
can give information about the mechanism that gener¬ 
ates a global current fluctuation. In general, rather than 
creating a global current fluctuation by a homogeneous 
contribution of particle jumps throughout the system, 
larger local currents are produced where the mobility is 
larger. For example, in the fluid regime of the ZRP this 
happens where more particles are available; in contrast, 
for the simple symmetric exclusion process we would ex¬ 
pect to see larger contributions to the global current for 
intermediate densities. It would be worthwhile to extend 
this picture to cases where there are dynamical phase 
transitions [38-40] leading to long-term accumulation of 
particles within the lattice. Further open questions relate 
to systems with non-diagonal diffusivity and mobility ma¬ 
trices (for the triangular lattice this can be achieved by 
setting p 2 P 3 )j as well as more general anisotropy with 
different physical processes in each direction. Experimen¬ 
tal tests of fluctuation relations in such situations would 


also be very enriching. 
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Appendix A: Derivation of the anisotropic 
fluctuation relation 

In this section, for completeness, we show a deriva¬ 
tion of the AFR [9] illustrating explicit details for the 
minimization of the action functional for current fluctua¬ 
tions. This fluctuation relation generalizes the IFR [8] to 
anisotropic systems, and was recently derived in [9] under 
the hypotheses that the system satisfies: i) reversibility 
and local detailed balance, ii) time-invariance of the ODP 
and OCP, and in) space-invariance of the OCP. We take 
as a framework the MFT to study systems satisfying the 
continuity equation 

^-V.i(r,() = 0, (Al) 

where p{r, t) and j{r, t) are the local particle density and 
local current respectively. Here we consider for the space 
variable, the d-dimensional unit interval H = [0,1]“^ which 
leads to the LDP for fluctuations of the global current 
as stated in Eq. (4). As in the main text, we consider a 
diffusive system in contact with two particle reservoirs 
with densities obeying pi > pr in the a:-direction, and 
periodic boundary conditions in every other direction. 

According to the MET, to compute the macroscopic 
RF we have to minimize [10, 41] 

e(J) = min - [ dr [ dr£ (r, r, p, Vp), (A2) 

Pd t Jq Jq 

with Lagrangian 

V7 ^ Uir,T)+D\'pf^{j{r,T) + D\'p) 

(A3) 

Here, we have that the local current is modelled by a 
deterministic and a stochastic term. The deterministic 
part, relates the current to the density via Tick’s law 
with diffusivity D(p) given by the diagonal matrix with 
elements Dk{p) = A^gip). Furthermore, the stochas¬ 
tic term corresponds to white noise with covari¬ 

ance L~‘^a{p)6{r' — r)6{t' — t). Here the mobility coeffi¬ 
cient is given by the diagonal matrix a{p) with elements 
crfc(p) = Efe(p)“^ = A^^/(p). Note that we have assumed 
the diffusivity and mobility matrices can be factorized as 
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a matrix of constant coefficients times a function of the 
density. The physical meaning of such a factorization is 
that particles diffuse at different rates in different direc¬ 
tions but obey a single type of process. In particular, if 
the constant matrices A and are both the identity, 
we have the isotropic dynamics for which the original IFR 
was derived. 

Since minimizing (A2) is still a very general problem, 
we now use hypotheses ii) time-invariant ODP and OCP, 
and iii) space-invariant OCP. Thus, the optimization 
problem is reduced to 

e(J) = mini / dr {J + D\7pf T. {J + DVp). (A4) 
p 4 


In contrast to [8, 9], we here explicitly solve the Euler- 
Lagrange equation 


dp 


E 


k=l 


d 

dxk 


dC 
. 9pi^^ 


= 0 , 


(AS) 


where we denote the space-variables in d-dimensions by 
Xk and p^J = d^'^^pjdx^ with k G {l,...,d}. Following 
this procedure one can compute that 


QH ^ [jk + A'fcpiV) PxIdpDk yl + ‘^JkDkpi} + Dl {pi^^ j dpak 


f)n 


dp 2au 


4^^ 


d \ dL \ dpDk j^p^^ldpDk Dkpx^dpak 

9xk I I cTfe 2(Tfc 2 ct? 


k [Sp 

Then, substituting these two expressions in (AS), some simplification leads to the differential equation 


E- 


2Dk [p^x^y dpDk + 2DIp^^} [dI (piV) - dpGu 


dCTk 




= 0 . 


(A6) 

(A7) 


(A8) 


Moreover, notice that periodic boundary conditions imply 
for the ODP that p^x} = 0 in all directions except for the 
one with open boundary conditions {xi), and allow us 

to replace 2Dlp^xi! by D\dp{p^xk)^■ Indeed, analogously 
to [42] in one dimension, we can now integrate Eq. (A8) 
with respect to the space variable. This yields the non¬ 
linear differential equation 


E 


k^l 



d 


E 


•^k 

dak 


+ C, 


(A9) 


where (7 is a constant of integration related to the bound¬ 
ary conditions. This way, to find the ODP that min¬ 
imizes (A4) we have to solve (A9), from which it can 
already be seen that for global current fluctuations lying 
on ellipses (constant first term on the RHS), the ODP 
will remain invariant. 

As a final step, note that Eq. (A9) can be written in a 
more compact way as 

{DWpf E (DVp) = J^EJ -h C. (AlO) 


It is easy to see that taking the difference between the 
RFs (A4) of two global current fluctuations for which 
the RHS of Eq. (AlO) has the same value, results in the 


relation 

e (J) - e (J') = l f dr (DVpf E (J' - J), (All) 
d r2 

where p is now the ODP. Here, corresponding to hypothe¬ 
ses ii) and iii) the OCP, J, can be taken out of the 
integral. Furthermore, from assumption i) it follows that 
the remaining integral in (All) is constant and we define 

E=l [ dr{D\/p)'^E. (A12) 

2 Jn 

This leads to the anisotropic version of Eq. (2), or in 
terms of the RF, 

e( J) - e( J') = E-{J' -J) (A13) 

for global currents satisfying 

J'^AJ = J'^AJ'. (A14) 

Here the density dependence, /(p), of the mobility ma¬ 
trix has cancelled out, and as expected this equation 
reduces to (3) for isotropic systems which satisfy the IFR. 
Furthermore, this symmetry also implies that the SCGF 
satisfies 


e(A) = e(A') 


(AI5) 
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for ellipses centred around the field E 

(A - Ef {\-E)= (A' - Ef A-i (A' - E) . 

(A16) 

In order to test this relation explicitly, one can calculate 


the ODP and the RF of the global current which has 
been done in 2-d for the KMP process [43] and the ZRP 
with interacting and non-interacting particles [9] (see also 
Appendix C). 


Appendix B: Quantum Hamiltonians for square and triangular lattices 


In this section we use the ladder operators (12) to write explicitly the Hamiltonians of the ZRP on the square 
and triangular lattices shown in Fig. 1. These Hamiltonians are equivalent to the stochastic generators with the 
corresponding geometry. At the boundaries, we use generic injection and extraction rates as shown in Fig. 7, whereas 
hopping rates for bulk sites are taken as shown in the insets of Fig. 1. For the square lattice we have 

L , 

- i7s = XI1“ (Sti - 1) + 7 (“Zi “ ^ (“jti - l) + - dj,L) 


L-1 

i=l 

i=i y 


(Bl) 


where assuming periodic boundary conditions in the j/-direction means that we identify j = L -I- 1 with j = I. Then, to 
measure current fluctuations, the Hamiltonian is modified by multiplying the terms corresponding to the bonds where 
we count particle jumps by the factors . Taking this into account, the modified Hamiltonian for the square lattice 
measuring current fluctuations globally and in the y-direction of region V (see Fig. 4) is given by 


e " - I 


- Hs = X (Xi® - 1 ) + 7 ^ (X 

L-1 


■ P ( l -. 


- 


— dj 


'J,L , 


2=1 

L 


-h 


XPy -f qy 


(X 


n _ (j 




(B2) 


Here, I(i, V) is the indicator function for the sites in V. 

For the triangular geometry the stochastic generator is given by 


E I O 

idr = X I X (Xi - 1) + 7fc {a~i - L.i) + 4 (aX - l) + 4 (a~L “ Lt 

L 

+ X^’i (ShX+i.i “ L.*) + 91 (XL”+I,i - dj+i,i) 


2=1 

M-1 


+ X [p^ (L.2iX2^-tl “ L.2i) + 92 (X2iS-.2*+l - dj,2i+l) 

2 = 1 

+ P3 (Qj,2iX+1.2i-|-l ~ L.2i) + 93 + ~ L + l,2i+l)] 


M 


X X (s+i.2*-iaX “ L-ti.2i-i) + 92 (aXi.2z-iaj.2i “ L. 2 O 


i=l 


+ P3 (aj.2*-lX2» - L.2i-l) + 93 (aXi-lL.2i - L.2*)] k 


(B3) 
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FIG. 7. Hopping rates for boundary sites of (a) square lattice and (b) triangular lattice. Input rates indicated in gray and 
output rates in black. 


whereas for the modified Hamiltonian counting 

L r 3 

— 

i = l U=2 

L 






1 . I cs 

= ^ ^ [ofe (a+ie”^'" - l) + 7fe (a“ie^'' - dj^i) + Sk + /3fe (ajLe~^^ - dj^L 

i=i U =2 

L 

2^1 
M-1 

[P^ («72*ajt2*+ie"^" - dj,2i) + 92 (ay2*a72i+ie^" “ dj^2i+i) 


2 = 1 


■P3 


M 


(®j,2i®^+l,2i+l® ^ dj^2i) + 93 (ay2i®j + 1.2i+l® ^ '^i + 1.2i+l)] 

'y2i^ ^ ~ '^i + 1.2i-l) + 92 ('2'j + l,2i-l®j,2i® ^ ~ dj^2i) 

(a 72 *-iS^ 2 ie”^" - di. 2 *-i) + 93 (ay 2 i-ia 72 ie^" - i^i. 2 *)] j- 


' [ 2*2 ('lj + l,2i-l®j,2i 

i=l 

+ P3 


(B4) 


Here we again assume periodic boundary conditions in the y-direction, and without loss of generality, an even number 
of sites L = 2M. 


Appendix C: Macroscopic RF and ODP 

In this appendix we show in detail how to calculate, 
according to the macroscopic fluctuation theory, the ODP 
and the RF for the 2-d ZRP. We follow [43], where a 
similar calculation was done for the 2-d KMP process. 
Again, we consider open boundary conditions in the x- 
direction and periodic in the y-direction. The left- and 
right-reservoir densities, pi and pr respectively, satisfy 
the inequality pi > Pr which indicates that the NESS 
has a mean current profile in the rightwards direction. 


We here assume the same hypotheses used to derive the 
AFR above (including space-homogeneous OCPs), solving 
Eq. (A4) for the RF and Eq. (A9) for the ODP. 


Firstly, note that the general diffusion and mobility 
(diagonal) matrices for the ZRP are given by D{p) = 
Az'{p) and a{p) = A-^z{p) (i.e., f{p) = z{p) and g{p) = 
z'ip)) where z'{p) = dz{p)/dp and the components are 
Afc = A^^ = pk [44, 45]. To compute the ODP, we 
substitute these in Eq. (A9) which leads to the non-linear 
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partial differential equation 


h <p) 


dp 

dxk 


V-^ + 4C. 

^iP^kZ(p) 


(Cl) 


Here, we denoted x by xi and y by X 2 . Note that the 
fugacity z{p) and its derivative z'{p) take different func¬ 
tional forms according to the type of interaction term Wn- 
Additionally, due to the cylindrical symmetry of the space 
we can assume that the profiles are flat in the direction of 
the periodic boundary conditions. Since we assume open 
boundary conditions in the cc-direction only, the k = 2 
term on the LHS of Eq. (Cl) vanishes leading to 


Px- 


Ap? 

z{p) 


dx 


2 

E 

k 


^^Px^z{p) 


4C. 


(C2) 


This equation is transformed into a differential equation 
for the fugacity, which depending on the type of interac¬ 
tion chosen, can be mapped to the density to obtain the 
optimal profile. 


with the negative sign corresponding to the physical solu¬ 
tion. 

In the non-monotonous regime the optimal profile has 
a maximum, z* = z(x*); the transition to this regime 
appears when the RHS of (C4) vanishes for the first time 
(i.e., when a;* = 0 and z* = zi). Moreover, from (C4) we 
identify b = —ajzi and from (C5) we get that the change 
of regime appears for currents 

— + — = 4:ZlPx {zi - Zr) . (C7) 

Px Py 

Clearly, these currents lie on an ellipse and will all have 
the same ODP. 

In the non-monotonous regime, we separate the solution 
of Eq. (C3) into two branches: one to the LHS of x* , and 
one to the RHS. Due to the non-linearity of this equation 
the derivative of the profile must be positive when x < x*, 
and negative when x > x* . Additionally, since z* is 
constant for any current on a fixed ellipse, we can use 
it to replace b by —ajz*. This way, we write the RHS 
of (C3) as a{\ — z/z*) finding that the ODP has fugacities 


1. Optimal density profile 


In the ZRP, the relation between the diffusivity and 
mobility matrices allows us to transform Eq. (A9) into a 
relation for the fugacity which can be solved independently 
of the interaction w^. The interaction plays a role when we 
relate the fugacity to the density, e.g., for non-interacting 
particles z{p) = p. Eor now, we eliminate the explicit 
dependence on the density and write Eq. (C2) in terms 
of the fugacity as 


^ , 4Cz 

dx) pI P^Py p^ 


(C3) 


The non-linearity of this equation requires us to con¬ 
sider two cases. The first one corresponds to a (de¬ 
creasing) monotonic ODP where the largest density is 
at the left reservoir, and the second, corresponds to a 
non-monotonous profile with a maximum fugacity, z*, at 
a distance x* from the left particle reservoir. 

Firstly, for the monotonous regime, we have to solve 

— = -\/a + bz, (C4) 

dx 

where a = J'llpl. + JyjpxPy and b = 4C/px- This leads 
to the solution 

. _ 5a;2 

z{x) = zi-x^Ja + bzl + (C5) 

where to satisfy the open boundary conditions, z(0) = z; 
and z(I) = Zri the constant of integration is determined 

by 


(C6) 


z{x) 




x<x* 
x > X* 


(C8) 

Here, the maximum fugacity z* and its position x* are 
determined self-consistently resulting in 


2 : = 


< (^Zi Zr \/n + 4:ZlZr'^ 

4(a- A2) 

, - A 2 {2zi + Va + 4:ZiZr) 
2{a-Al) 


(C9) 


where A^ = zi — Zr- 


2. Global current rate function 


In addition to the ODP, we can calculate exactly the RE 
of the ZRP. This means solving Eq. (A4) constrained by 
Eq. (AlO) which, by the symmetry of our system reduces 
to calculating the integral 



(^Jx + z'{p)px j2 

4pxz{p) 4pyz(p) 


dx, 


(CIO) 

with the minimizing constraint (C2). Similarly to the 
procedure above, we use Eq. (C3) to work in terms of 
the fugacity (instead of the density), and substitute in 
Eq. (CIO) to calculate the RE. Note that we still have 
to account for the change of regime for currents larger 
than the threshold given in (C7). This means that, in the 


6 = 4 (z; -I- Zr- ± + 4z;Zr.) 
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monotonous scenario, the RF is obtained by integrating and (C12) is given by 


e(J) 


+ j.) 


(Cll) 


^(j) = Y (f^) T ("^0 + fei + Va + 

— 2-\/osinh“^ ~ 2-\/asinh“^ 

" (C13) 

for the monotonous regime, and 


while, for the non-monotonous regime, it becomes 



Px^ (l 2z* ) j 

PxO- (l ~ 2IV) _ Jx 


dz 


dz 


(C12) 


where the constants C and z* are determined by (C6) 
and (C9). The exact solution of the integrals (Cll) 


KJ) = 


- In 


2 


+ Px'/a 


(C14) 

for the non-monotonous regime. For each specific ZRP 
interaction, the corresponding relation between density 
and fugacity can be used to obtain the RF in terms of the 
reservoir densities. In particular, for the case = n we 
can simply replace z by p, whereas for Wn = w we replace 
z by pj{p -f 1). Finally, to compare with the microscopic 
approach we have to relate the boundary fugacities to the 
boundary rates as mentioned in the main text. 
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